0. Why Parallelise a Simulation?

1. Capability — the model does not fit in one machine's memory (a 10¹⁰-cell mesh, a 10⁸-agent population).
2. Capacity / turnaround — the run must finish before the decision is needed (a weather forecast for tomorrow is useless if it takes 30 hours).
3. Statistical demand — the 1/√n law of Unit I means accurate stochastic results need thousands of replications.
4. Geographic distribution — federated simulations (HLA) join simulators owned by different organisations that cannot be co-located.
5. Hardware reality — single-core clock speed stopped increasing around 2005; all further performance is parallel.

0.1 The limits: Amdahl and Gustafson

Let f be the serial fraction of the work and p the number of processors.

With only 5% serial code, no machine on earth can exceed a 20× speedup on a fixed problem — that is Amdahl's warning. Gustafson's answer is that in practice we do not keep the problem fixed: given more processors we simulate a finer mesh or a bigger population, and the serial fraction shrinks relative to the growing parallel work. Hence the two notions of scaling in Section 6.

1. Partitioning the Data

This is the dominant strategy in scientific simulation because it scales with problem size, whereas the number of distinct algorithmic stages is fixed.

1.1 Geometric decompositions of a mesh

1.2 Halo (ghost) regions

A stencil computation needs the values of neighbouring cells that live on another processor. Each partition therefore keeps a halo — a read-only copy of the boundary layer of its neighbours — that is refreshed by a halo exchange once per step. Halo width equals the stencil radius; a wider stencil or a higher-order method costs proportionally more communication.

for each time step:
    start_halo_exchange()          # non-blocking Isend/Irecv
    compute_interior()             # overlap communication with computation
    wait_for_halo()
    compute_boundary()
    swap(old, new)

1.3 Partitioning other data structures

• Particles — partition by spatial region (particles migrate between processors as they move, so ownership changes) or by index (simple but destroys locality).
• Agents — partition the environment; agents crossing a boundary are serialised and sent to the new owner.
• Graphs — partition vertices to minimise edge cut (Section 7).
• Replications — the embarrassingly parallel case: each processor runs a complete, independent replication with its own random substream. Near-linear speedup, no communication. Always do this first if the goal is statistical accuracy rather than a single big run.

2. Partitioning the Algorithms

2.1 Forms of algorithmic parallelism

1. Pipelining. Stages arranged in sequence, each working on a different time step or dataset. Throughput rises by the number of stages; latency does not. Efficiency requires balanced stage times, and there is a fill/drain overhead.
2. Component (model) parallelism. In a coupled climate model, the atmosphere, ocean, land and ice components run on separate processor groups and exchange fluxes at coupling intervals.
3. Task parallelism with a task graph. Work is expressed as a DAG of tasks with dependencies; a runtime (Charm++, StarPU, Dask, Ray) schedules ready tasks onto free workers. This handles irregular workloads far better than static partitioning.
4. Parameter sweep / replication parallelism. Independent runs across the scenario space — a task-parallel pattern with zero communication.
5. Solver-internal parallelism. Parallel linear algebra, parallel FFT, parallel sorting inside a single time step.
6. Speculative execution. Compute a branch before knowing whether it is needed; the basis of optimistic PDES (Section 3.4).

2.2 Data vs. algorithm partitioning

Real large-scale codes use both: components are assigned to processor groups (algorithmic), and each group decomposes its own domain (data). That is the hybrid partitioning of Section 8.

3. Handling Inter-Partition Dependencies

Cutting a model creates dependencies across the cut. Managing them correctly — without destroying performance — is the central technical content of parallel simulation.

3.1 Kinds of dependency

• Spatial (boundary) dependency — a cell needs neighbour values. Handled by halo exchange.
• Temporal (causal) dependency — an event on partition A at time t may cause an event on B at t+δ. Handled by synchronisation protocols.
• Global (reduction) dependency — a shared quantity (total energy, a global time step, a convergence test) needs all partitions. Handled by collective operations, and it is the main source of the serial fraction in Amdahl's law.
• Resource dependency — entities on different partitions compete for one shared resource; needs a distributed protocol or resource ownership.

3.2 Synchronous (time-stepped) execution

The easy case. All partitions execute step n, exchange haloes, hit a barrier , and proceed to step n+1. Correctness is automatic; the cost is that every barrier runs at the speed of the slowest processor, so load imbalance and OS jitter accumulate.

3.3 Conservative synchronisation

An LP processes an event only when it can prove that no earlier event can still arrive.

1. Each incoming link carries a clock: the timestamp of the last message received on it.
2. The LP may safely process any event with timestamp less than the minimum over all input links.
3. If a link is silent the LP blocks → potential deadlock.
4. Chandy–Misra–Bryant null messages break the deadlock: an LP that will send nothing before time t sends a null message with timestamp t ("nothing from me until t").
5. Performance depends critically on lookahead — the guaranteed minimum delay before an LP can affect another. Zero lookahead means no parallelism at all.

3.4 Optimistic synchronisation (Time Warp)

1. Each LP processes events as fast as it can, assuming no straggler will arrive.
2. State is saved (checkpointed) periodically so that it can be restored.
3. When a straggler with timestamp t_s arrives, the LP rolls back to the last state before t_s.
4. Messages sent in error are cancelled by anti-messages, which may cause cascading rollbacks in other LPs.
5. Global Virtual Time (GVT) — the minimum timestamp of any unprocessed event or message in flight — is computed periodically. Nothing before GVT can ever be rolled back, so memory for older checkpoints is reclaimed (fossil collection) and irrevocable actions such as I/O are committed only up to GVT.

3.5 Reducing dependency cost

• Overlap communication with computation using non-blocking calls.
• Aggregate many small messages into one large one — latency, not bandwidth, usually dominates.
• Replicate cheap read-only data instead of fetching it.
• Relax where the model permits: an asynchronous or slightly stale value may be acceptable, and the error should then be quantified.
• Increase lookahead by exploiting known minimum delays in the physical system.

4. Dynamic Partitioning and Load Balancing

4.1 Why the load moves

• Adaptive mesh refinement adds cells where a shock or flame front travels.
• Particles or agents cluster: a crowd gathers, galaxies collapse, traffic congests.
• Physics changes cost per cell: chemistry activates only in the reacting region.
• Heterogeneous or shared hardware: some nodes are slower, or shared with other jobs.
• Faults reduce the resource pool mid-run.

4.2 The four questions of any load-balancing scheme

1. Measure — what is the load metric? Cells, particles, events processed, or measured wall-clock time per rank (usually the most honest).
2. Decide when — every k steps, or when imbalance λ = max/mean exceeds a threshold. Balance the benefit of rebalancing against its cost.
3. Decide how — diffusive (shift work to lighter neighbours; small data movement, slow convergence) or global repartitioning (recompute the whole partition; better quality, more movement).
4. Migrate — serialise the state, transfer, rebuild indices and neighbour lists, resume. Migration must preserve correctness of in-flight messages.

4.3 Techniques

• Over-decomposition + work stealing. Create many more partitions than processors and let idle workers steal. Used by Charm++, Cilk, task runtimes; adapts automatically and tolerates heterogeneous hardware.
• Master–worker (bag of tasks). Excellent for parameter sweeps and replications; the master becomes a bottleneck at very large scale, fixed by hierarchical masters.
• Space-filling curve repartitioning. Reorder entities along a Hilbert curve and cut into equal-work segments — cheap, incremental, locality-preserving; the standard in AMR and N-body codes.
• Repartitioning graph tools — ParMETIS, Zoltan, Scotch; these explicitly optimise the sum of edge cut and migration cost.
• Predictive balancing — use the previous steps' measurements to forecast the next; effective because simulation workloads evolve smoothly.

5. Communication Patterns in Partitioned Systems

On a modern cluster α ≈ 1–5 µs while 1/β is tens of GB/s. Sending 1000 messages of 8 bytes costs a thousand latencies; sending one message of 8000 bytes costs one. This is why message aggregation is the first optimisation to try.

5.1 The standard patterns

5.2 Interest management

In distributed simulation with many entities, sending every update to everyone is O(n² ) and impossible beyond a few thousand entities. Interest management (HLA Data Distribution Management) filters updates so that a federate receives only what it cares about, using routing spaces: the world is divided into regions or grid cells, each entity publishes to the cells it occupies and subscribes to the cells it can perceive. This turns a global broadcast into localised multicast.

5.3 Practical rules

1. Prefer non-blocking communication and overlap it with interior computation.
2. Aggregate small messages; pack contiguous buffers rather than sending strided data.
3. Remove unnecessary collectives — a global reduction every step is often needed only every tenth step.
4. Keep the communication topology matched to the network topology where possible (rank reordering).
5. Measure with a profiler (Score-P, TAU, mpiP, Nsight) before optimising; intuition about where the time goes is usually wrong.

6. Scalability Challenges in Partitioned Systems

6.1 The eight barriers to scalability

1. Serial fraction — initialisation, I/O, global decisions (Amdahl).
2. Load imbalance — the barrier runs at the slowest rank; even 5% imbalance caps efficiency near 95% and it compounds every step.
3. Communication growth — in strong scaling the partition shrinks, so the surface-to-volume ratio worsens and communication eventually dominates.
4. Synchronisation and jitter — the barrier amplifies rare OS or network hiccups across all ranks.
5. Collective operations — O(log p) at best, and allreduce latency becomes the floor of the time step at very large p.
6. Memory per node — replicated global structures (a full agent directory, a full mesh copy) do not shrink with p and eventually exhaust memory.
7. I/O bottleneck — thousands of ranks writing checkpoints to one file system; solved by parallel I/O (MPI-IO, HDF5, ADIOS) and in-situ analysis (Unit V).
8. Fault probability — with more components, mean time between failures falls; at extreme scale the run will be interrupted (Section 9).

6.2 Diagnosing scalability

• Plot speedup and efficiency against p on log axes; the departure point identifies the limit.
• Instrument the phases: computation, communication, idle/barrier wait, I/O. Growing idle time means imbalance; growing communication means the partition is too small or too badly shaped.
• Use a roofline or a simple T(p) = T_comp/p + T_comm(p) + T_serial model to predict where the curve will bend, before buying time on a larger machine.

7. Partitioning in Graph-Based Systems

7.1 Why graphs are the hard case

A mesh has geometry, so a coordinate cut works. A social, citation or web graph has no useful geometry, has a power-law degree distribution, and has a small diameter — so every balanced cut separates many edges. This is why distributed graph processing is dominated by communication.

7.2 Methods

• Multilevel partitioning (METIS/ParMETIS, Scotch, KaHIP) — the standard and the one to name in an exam: coarsen the graph by matching and merging vertices, partition the tiny coarse graph well, then uncoarsen and refine at each level with Kernighan–Lin / Fiduccia–Mattheyses local swaps.
• Spectral partitioning — use the Fiedler vector (eigenvector of the second smallest eigenvalue of the Laplacian) to order and split vertices. High quality, expensive.
• Geometric — recursive coordinate bisection, space-filling curves; applicable only when vertices have coordinates.
• Streaming / one-pass heuristics (LDG, FENNEL) — assign each vertex as it arrives; necessary when the graph does not fit in memory.
• Label propagation and community detection — exploit natural cluster structure; effective on social graphs.

7.3 Vertex-cut versus edge-cut

For power-law graphs, cutting vertices beats cutting edges. In the vertex-cut model (PowerGraph/GraphX) a high-degree hub is replicated across the machines that hold its edges, and its state is reconciled by a small gather–apply–scatter step. Because a hub with a million edges cannot live on one machine without wrecking the balance, this reformulation was the key advance for real-world graph systems.

7.4 Programming models for graph simulation

• BSP / Pregel — "think like a vertex": supersteps of local computation, message exchange, then a global barrier. Simple and deterministic; barrier-bound.
• GAS (gather–apply–scatter) — PowerGraph's model, designed around vertex cuts.
• Asynchronous — no barrier; faster convergence for iterative algorithms but non-deterministic, which conflicts with simulation reproducibility.

8. Hybrid Partitioning Approaches

8.1 Forms of hybridisation

1. MPI + OpenMP + CUDA. MPI between nodes (distributed memory), threads within a node (shared memory), and GPU kernels for the inner loops. Fewer, larger MPI partitions means less surface area and fewer messages.
2. Data + task hybrid. Domain decomposition for the mesh, plus a task runtime for the irregular parts (chemistry, particle physics, adaptivity).
3. Component + domain hybrid. Coupled multiphysics: each component owns a processor group (algorithmic partitioning) and decomposes its own domain (data partitioning).
4. Multi-constraint / multi-objective partitioning. Balance two loads at once — e.g. cell count and particle count in a particle-in-cell code — which single weighted partitioning cannot do.
5. Static + dynamic hybrid. A good static partition at start-up, with cheap diffusive corrections during the run.
6. Mixed synchronisation. Conservative synchronisation between well-separated clusters with good lookahead, optimistic within a tightly coupled cluster.

8.2 Why hybrid usually wins

• It matches the memory hierarchy: intra-node communication through shared memory is far cheaper than through the network.
• It reduces the number of partitions, which improves the surface-to-volume ratio and shrinks halo volume.
• It lets each part of a heterogeneous model use the decomposition that suits it.
• It exploits accelerators for the regular inner kernels while leaving irregular control flow on the CPU.

9. Fault Tolerance in Partitioned Systems

At 10⁵ nodes, even a per-node MTBF of ten years gives a system failure roughly every hour. A simulation that runs for a day must therefore be able to survive failures.

9.1 Failure model

• Fail-stop — a node dies and stops. The common assumption.
• Silent data corruption — a bit flips and the run continues with wrong numbers. The dangerous case, detected by residual checks, invariants or replication.
• Network partition — groups can no longer communicate; relevant to distributed (rather than tightly coupled) simulation.
• Straggler / slow node — not a failure but has the same effect at a barrier.

9.2 Checkpoint and restart

Techniques that reduce C: incremental checkpoints (save only changed pages), multilevel checkpointing (node-local SSD → partner node → parallel file system, as in SCR/FTI), asynchronous checkpointing (write in the background), and compression. In an application-level checkpoint the model writes only the physical state it actually needs, which is usually far smaller than a system-level memory image — and it is portable across machines.

9.3 Beyond checkpointing

• Replication — run duplicate partitions; expensive (2× or more) but detects silent corruption by comparison and eliminates restart delay.
• Algorithm-based fault tolerance (ABFT) — encode checksums into the data structures (rows/columns of a matrix) so that a lost or corrupted block can be reconstructed arithmetically.
• Natural fault tolerance — iterative and stochastic simulations can sometimes absorb the loss of a partition with a quantified error, rather than restarting.
• Message logging — log messages and replay them so that only the failed process rolls back, not the whole job; the natural fit for PDES, where rollback machinery (Time Warp) already exists.
• Resilient runtimes — ULFM MPI, Charm++ and workflow engines that detect failure, shrink or replace the process set and rebalance.
• Task-level retry — in a replication or parameter sweep, a failed run is simply re-queued. This is why embarrassingly parallel work is also the most robust.

10. Partitioning for Emerging Architectures

The right partition depends on the machine. As hardware diversifies, partitioning must become architecture-aware.

10.1 Cross-cutting trends

1. Data movement, not arithmetic, is the cost. Moving a double across a node can cost more energy than a hundred floating-point operations, so partitioning is increasingly an energy-optimisation problem.
2. Asynchrony over barriers. Task-based, dependency-driven runtimes tolerate heterogeneity and jitter better than bulk-synchronous code.
3. Performance portability. Kokkos, RAJA, SYCL and OpenMP target-offload let one partitioned code run on CPU and multiple GPU vendors.
4. In-situ analysis. Because I/O does not scale, analysis and visualisation move into the simulation itself — the bridge to Unit V.
5. Surrogates and ML acceleration. Learned models replace expensive kernels, and the partitioning must then balance a mix of physics and inference work.

11. A Practical Conversion Recipe

If an examiner (or a project) asks “how would you convert this sequential simulation to a parallel one?”, answer with this sequence:

1. Profile first. Find the hot loops and the memory footprint; do not parallelise what does not matter.
2. Try replication parallelism first. If the goal is statistical accuracy, run independent replications — near-linear speedup for almost no work.
3. Choose the decomposition. Data if the state is large and homogeneous; algorithmic if the model is a pipeline of distinct components; hybrid at scale.
4. Identify dependencies. Spatial (halo), temporal (causality), global (reductions), resource contention.
5. Choose a synchronisation scheme. Barrier per step for time-stepped models; conservative or optimistic for event-driven ones, chosen by available lookahead.
6. Design the communication. Non-blocking, aggregated, overlapped; remove unnecessary collectives.
7. Handle randomness. Independent substreams per partition, and ensure the result is independent of the processor count.
8. Verify against the sequential run. Same seed, same answer (or the same statistics with a documented reason for bitwise differences).
9. Measure strong and weak scaling, and attribute lost efficiency to imbalance, communication or serial work.
10. Add resilience — checkpoints at the Young/Daly interval, and a restart path that has actually been tested.

12. Unit Summary

• Parallelisation is a decomposition problem: balance the work, minimise the cut.
• Amdahl bounds strong scaling by 1/f; Gustafson explains why weak scaling still works.
• Data (domain) partitioning splits the state and runs the same code everywhere; compactness matters because computation scales with volume and communication with surface.
• Algorithm (functional) partitioning splits stages or components; limited scalability, but essential for coupled multiphysics and pipelines.
• Dependencies are spatial (halo exchange), temporal (causality constraint), global (collectives) and resource-based. PDES resolves temporal dependencies conservatively (lookahead, null messages) or optimistically (rollback, anti-messages, GVT).
• Dynamic partitioning answers four questions: measure, when, how, migrate. Over-decomposition with work stealing and space-filling-curve repartitioning are the practical techniques.
• Communication cost is α + nβ; aggregate messages, overlap with computation, avoid all-to-all and unnecessary reductions, and use interest management in distributed simulation.
• Scalability is limited by serial fraction, imbalance, communication growth, jitter, collectives, memory, I/O and faults.
• Graph partitioning is NP-hard; multilevel methods (METIS/ParMETIS) are standard, and vertex cuts handle power-law graphs better than edge cuts.
• Hybrid partitioning matches the hardware hierarchy (MPI+OpenMP+GPU) and mixes strategies and synchronisation schemes.
• Fault tolerance at scale: checkpoint/restart at the √(2CM) interval, multilevel and incremental checkpoints, ABFT, replication, message logging and resilient runtimes.
• Emerging architectures shift the objective from balancing FLOPs to minimising data movement and energy.

12.1 Key terms

SPMD/MPMD · domain decomposition · halo/ghost cells · surface-to-volume ratio · Amdahl and Gustafson · strong and weak scaling · logical process · causality constraint · straggler · lookahead · null message · Time Warp · anti-message · GVT · fossil collection · work stealing · over-decomposition · space-filling curve · edge cut · vertex cut · multilevel partitioning · BSP superstep · interest management · allreduce · checkpoint interval · ABFT · ULFM.

12.2 Practice questions

Short answer (2–3 marks each)

1. State Amdahl's law and give the maximum speedup when the serial fraction is 10%.
2. Differentiate strong scaling from weak scaling.
3. What is a halo region and why is it needed?
4. Define lookahead and explain its role in conservative synchronisation.
5. What are anti-messages in Time Warp?
6. Why is a vertex cut preferred to an edge cut for power-law graphs?
7. Write the optimal checkpoint-interval formula and define its symbols.

Medium answer (5 marks each)

1. Compare data partitioning and algorithm partitioning under at least five criteria.
2. Explain the surface-to-volume principle and use it to compare 1-D, 2-D and 3-D decompositions of a grid.
3. Describe the four decisions involved in dynamic load balancing and two techniques used in practice.
4. List the standard communication patterns with their cost in p and one simulation use for each.
5. Explain multilevel graph partitioning (coarsen – partition – uncoarsen/refine).
6. Describe three fault-tolerance techniques for large parallel simulations other than plain checkpointing.

Long answer (10 marks each)

1. Explain the causality problem in parallel discrete-event simulation and compare conservative and optimistic synchronisation in detail, with examples of when each is preferred.
2. Discuss the scalability challenges of partitioned simulations, with an analytical example showing why strong scaling saturates, and describe how each challenge is mitigated.
3. Given a large agent-based epidemic simulation on a contact network, describe in full how you would convert it to a parallel and distributed implementation: partitioning, dependencies, synchronisation, communication, load balancing, randomness, verification and resilience.
4. Discuss partitioning strategies for emerging architectures (GPU, heterogeneous nodes, cloud, edge), explaining how the objective changes from balancing computation to minimising data movement.
5. Explain hybrid partitioning approaches and justify, with reasons and counter-arguments, when MPI+OpenMP+GPU is preferable to a flat MPI decomposition.

12.3 Further reading

• R. M. Fujimoto, Parallel and Distributed Simulation Systems, Wiley — the standard text for Sections 3 and 9.
• G. Karypis and V. Kumar, METIS / ParMETIS technical reports — for Section 7.
• A. Grama, A. Gupta, G. Karypis and V. Kumar, Introduction to Parallel Computing — for Sections 1, 2, 5 and 6.
• J. Banks (ed.), Handbook of Simulation — chapter on parallel and distributed simulation.
• J. Daly, “A higher order estimate of the optimum checkpoint interval”, Future Generation Computer Systems, 2006.